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Abstract 

This paper concerns the reconstruction of an anisotropic diffusion tensor 7 = ("fij)i<i,j<2 
from knowledge of internal functionals of the form jVui ■ Vitj with Ui for 1 < i < I solutions 
of the elliptic equation V • 7V?/.; = on a two dimensional bounded domain with appropriate 
boundary conditions. We show that for I ~ A and appropriately chosen boundary conditions, 
7 may uniquely and stably be reconstructed from such internal functionals, which appear in 
coupled-physics inverse problems involving the ultrasound modulation of electrical or optical 
coefficients. Explicit reconstruction procedures for the diffusion tensor are presented and 
implemented numerically. 

1 Introduction 

Coupled-physics modalities are being extensively studied in medical imaging as a means to com- 
bine high contrast with high resolution. Such imaging modalities typically couple a high-contrast 
low-resolution modality with a low-contrast high-resolution modality. In this context, Ultra- 
sound Modulated Optical Tomography (UMOT) or Ultrasound Modulated Electrical Impedance 
Tomography (UMEIT) aim to improve the low resolution in the reconstruction of diffusion (in 
OT) or conductivity (in EIT) coefficients by perturbing the medium with focused or delocalized 
ultrasound when making measurements. We assume here that the electric potential in EIT and 
the photon density in OT are modeled by the following elliptic model 

n 

-V • (jVu) = ~ ^2 di ilijdju) = 0, in X, u\ ax =9 on dX, (1) 
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where X is a subset in W 1 , where n will equal 2 in this paper, and where 7 is a symmetric 
positive definite tensor. In the case of OT, equation ([1]) is an approximation of a more accurate 
model that takes into account absorption effects by adding a zeroth order term a a u with a a > 
in the left-hand side of (pQ). Whether this additional term can be handled with the present 
approach will be the object of future research. We also assume that the ultrasound perturbation 
of the medium and of the corresponding boundary (current) measurements allow us to stably 
reconstruct the power density of two solutions u, v of ([1]) with prescribed boundary conditions, 
namely, the quantity H[u,v] := 7VU • over X. How to obtain power densities in practice 
has been addressed in e.g. p] using highly focused ultrasonic waves, and in e.g. p2J [7] in the 
context of synthetic focusing. 

The main objective of this paper is the reconstruction of the symmetric tensor 7 in ([I]) and 
in dimension n = 2 from knowledge of internal measurements {H[v,i, Uj]}Y) = i, where each Uj 
corresponds to a different, properly chosen boundary illumination g^. 

The isotropic case 7 = al n was analyzed in two dimensions in [10] and extended to three 
dimensions in [7] and to all dimensions n > 2 and more general types of measurements of the 
form o 2a \Vu\ 2 ,a £ R, (n - 2) a + 1 / in [16]. The case of internal current densities of the 
form |7Vu| arises in current density impedance imaging, see e.g. [T7] and [6] for a review and 
bibliography of recent results obtained in similar problems. Although the results in this paper 
also extend to general values of a, we restrict ourselves to the case a = ^ to simplify the 
presentation. Similar problems in dimensions n = 2 and n = 3 were also analyzed in a linearized 
context in [13]; see also the recent paper p3] on general linearized hybrid inverse problems. 

This paper extends the analyzes in [16] to the case of two-dimensional anisotropic diffusion 
tensors. We show that with 4 (or, in practice, in fact 3) properly chosen illuminations, the mea- 
surements {Hij}jj =1 allow for a unique and stable reconstruction of the full anisotropic tensor 
7. In particular, the internal functionals considered here do allow us to uniquely reconstruct 
the conformal structure (or normalized anisotropy tensor) 7 := (det 7) "2^ unlike the case of 
Dirichlet-to-Neumann boundary measurements as they appear in the Calderon problem, where 
7 can be reconstructed up to a push-forward by an arbitrary change of variables leaving each 
point on dX invariant [3]. 

More precisely, using the decomposition 7 = (det 7) 27 with det 7 = 1, the anisotropy tensor 7 
can be reconstructed in an algebraic and pointwise manner, after which the quantity (det 7) 2 can 
be obtained in two possible ways, either via inverting two consecutive gradient (or, after taking 
divergence, Poisson) equations, or by inverting a strongly coupled elliptic system followed by a 
gradient (or Poisson) equation. The reconstruction of (det 7) 2 is similar to the case where 7 is 
known up to multiplication by a scalar function. Since this problem has the same dimensionality 
as that of the reconstruction of an isotropic diffusion tensor (treated in [101 \7\ [16] ) , only m = n 
illuminations are necessary and the reconstruction can be done following the second step of the 
previously described approach. Although some of the techniques presented here generalize to 
higher dimensions, we restrict ourselves to the two-dimensional setting in this paper. 
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Finally, some numerical simulations confirm the theoretical predictions: both the isotropic 
and the anisotropic parts of the tensor can be stably reconstructed, with a robustness to noise 
that is much better for the former than the latter. 

Our main results are presented in section[2l Their derivation occupies sections [3][5l Numerical 
simulations are shown in O and concluding remarks offered in section [7J 

2 Statement of the main results 

Let 1 C I 2 be an open, bounded, simply connected domain. Borrowing notation from [2], for 
k > 1, a real 2x2 symmetric matrix 7 = {jij} 2 j = i belongs to M S K (2) if and only if it satisfies 
the uniform ellipticity condition 

K~ l \i\ 2 < -fijCe < k|£| 2 for every f G R 2 . (2) 

In the following, we consider the problem of reconstructing an anisotropic conductivity function 
7 G L°°(X,A4 S K0 (2)) where kq > 1 is fixed, from the knowledge of power-density measurements 
of the form 

Hij(x) = jVui • Vuj, 1 < i, j < m, 

where for 1 < i < m, the function m satisfies the conductivity equation 

-V ■ (jVui) = in X, Ui\ ax = 9i ondX. (3) 

The gfj's are prescribed illuminations. We denote by A the unique positive A^ s (2)-valued function 
that satisfies the pointwise relation A 2 (x) = 7(x). Clearly, A E L°° (X , Ai s ^—(2)) . 

We first change unknown functions by defining for 1 < i < m, the vector fields Si := AVui. 
The elliptic equation ([3]) thus reads 

V • (AS,) = 0. (4) 

Furthermore, denoting the "curl" operator in 2D by JV- where J = the fact that 

A~ l Si = VtXj implies that it is curl-free, that is, 

JV • (A" 1 Si) = 0. (5) 

The data become Hij = Si • Sj. We now decompose A into 

A = \A\^A, \A\:=detA and detl=l. (6) 

From the uniform bounds k 2 < \A\i < k§ , it follows immediatly that A G L°°(X,A4 S ). 
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Over any open, connected set fl Q X, where two solutions 51,5*2 satisfy the following 
positivity condition 

inf (det = inf det(S 1 ,S 2 ) > c > 0, H = {S t ■ 5,}|, =1 , (7) 

we can derive the first important relation 

V log |.A| = ^Vlog|tf| + (VH jl ■ ASi)A- l Sj = \H\-^(V(\H\^W l ) ■ AS^A^Sj. (8) 

We now orthonormalize the frame S = (Si,S 2 ) into a 50(2)-valued frame R = (Ri, R 2 ), via 
a transformation of the form Ri = UjSj (or, in matrix notation R = ST T , T := {tij}), where 
T is known from the data. For further use we denote T~ 1 = and define the vector fields 
{Vij}j =1 as 

Vy := Vfefcjt*'', l<i,j<2, V<*:=^(V tJ -V jt ). (9) 

This orthonormalization can be obtained by the Gram-Schmidt procedure or by setting T = 
for instance. Since R is >SO(2)-valued, we parameterize it with a S 1 -valued function 9 such that 
^ = [sine "cosfl ] ' ^ ne * s then able to derive the following second important equation by using 
(|4|) and geometric arguments: 

A 2 V9 + [A 2 , M] = A 2 V& - ^JN, (10) 

where N := log \H\ and Vf 2 are known from the data, and [A 2 , A\] := (A 2 -'V)Ai — (Ai-V)A 2 , 
where Aj denotes the j-th column of A. 

Reconstruction of the anisotropy A: We start by defining the set of admissible illumina- 
tions for some 7 E L°°(X, M S K (2)), by stating that a quadruple g = (g±, g 2 , gs, 54) belongs 
to Q-y if the following conditions are satisfied for some constant cq > 0: 

inf_min(det(Vui, Vu 2 ), det(Vu 3 , Vm 4 )) > c > 0, (11) 

xex 

Y s := -JVlog(det(Vui, Vu 2 )/det(Vu 3 , Vu 4 )) + for every x G X, (12) 

where u% solves ([3]) for 1 < % < 4. Condition (1121) . which is directly motivated by calcula- 
tions later, expresses the fact that the relative variations of the volumes det(Vtti, VU2) and 
det(Vn3,Vn4) differ at every point, which seems to be the condition that guarantees that our 
measurements are rich enough to "see" the anisotropy. 

Condition (|lip is rather easy to ensure by virtue of [U Theorem 4], which guarantees that 
(TTT1) holds as soon as both maps dX B x 1— > (gi(x), g 2 (x)) and dX B x 1— > (53 (x), 54 (x)) are 



4 



homeomorphisms onto their images. Based on a construction that uses Complex Geometrical 
Optics (CGO) solutions, we construct in the next lemma solutions that satisfy condition (|12p 
under some regularity assumption on 7. This in turn guarantees that C/ 7 is not empty when 7 
is smooth enough. 

Lemma 2.1. Let 7 G L°°(X,A4 s Ko (2)) for some kq > 1 be such that |-y| 2 G H 5+£ (X) for some 
e > and the function v : X — > C defined as 

A 3 1 1 — y u(x) = -(x) (13) 

7ii + 722 + 2|-y| 2 



is locally of class C over X. Then the set Qy defined by conditions (|lip and (| 12 j) is not empty 
and contains an open set of sufficiently smooth boundary conditions. 

Consider now (g 1 , g 2 , g 3 , 9a) G Qj, and let us denote gW = (51,92) and g( 2 ) = {g 3 ,gi). 

Define two orthonormal frames = [B,\ \R 2 ], i = 1,2, obtained after orthonormalization of 

gW = [5i|52], and = [S3IS4]. Taking the difference of equations (fTUj) for each system, we 
obtain the algebraic equation 

A 2 X S = Y" g , where X s := V(6 2 - 0i) - V$ 2) + V£ (1) , (14) 



and Yg is defined in (|12p . Both vector fields X s and Yg are uniquely determined by the data, 
see below. Since A is only described by two scalar parameters, equation (|14p together with the 
condition (|12p allow us to reconstruct these two parameters algebraically and in a pointwise 
manner. When orthonormalization uses the Gram-Schmidt procedure, X s and Yg satisfy the 
following boundedness and stability inequalities for some constants C\ , C 2 ■ 

max(\\X g \\ L °°( X ), \\Y\\ L oo( X }) < Ci\\H\\ w i,oo( X ), 
max(||X g - X' s \\ L <x,( X ), ||Y g - Yg\\ L <x>( X )) < C 2 \\H - H'\\ w i,°o( X ), 

where H = {Hij}fj =1 and H' = {Hlj}fj =l respectively come from g G C? 7 and g G Gy. We 
arrive at the following result: 

Theorem 2.2 (Anisotropy reconstruction in 2d with m = 4). For g G Q^, the measurements 
H = {Hij}fj =l uniquely determine the tensor A via the explicit algebraic equation (|14p . More- 
over, for 7,7' withg G Qy andg G Qy with the corresponding measurements Hij,H-j G W 1,0O (X) 
for 1 < i,j < 2, and in the case where orthonormalization is done using the Gram-Schmidt pro- 
cedure, the following stability statement holds: 

\\A-A!\\ Lao{kX ) < C\\H -H'\\ w i,oo, (16) 

for some constant C. 

Remark 2.3. In practice, we have observed numerically that m = 3 was enough to reconstruct 
7 if we chose g G £ 7 of the form {gi, g 2 , g 2 , g 3 ). 
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Reconstruction of (#, log or (u\,U2, \ A]" 1 ): Once the anisotropy A is known, the prob- 
lem of reconstructing \A\ now has the same dimensionality as that of reconstructing an isotropic 
conductivity. It requires only m = 2 internal functionals satisfying ([7]) in X. 

A first approach towards reconstructing \A\ consists in solving a gradient equation for the 
scalar quantities 9 and then log the right-hand sides of which are successively known. If 
9 and \A\ are known throughout the domain's boundary, one may take the divergence of said 
gradient equations instead and solve Poisson equations with Dirichlet boundary conditions. 

Such an approach provides Lipschitz-stable reconstructions as is summarized in 

Theorem 2.4 (Stability of \A\, first approach). Assume that A is either known or recon- 
structed as in theorem \2.2\ with the stability estimate (|16p . Then \A\ is uniquely determined 
by {Hij}i<ij<2 £ W 1,0 ° satisfying ([7]). Moreover, if two set H and H' jointly satisfy the pre- 
vious assumptions, the corresponding reconstructed coefficients \A\ and \A'\ satisfy the stability 
inequality 

II 1°§ I ^4 1 — log \A'\ || w 1 '°°{x) 

< C || H-H' || wl ,„. (17) 

A second approach consists in inserting the expression in equation (|8|) into the elliptic equa- 
tion ([3]) and deriving a strongly coupled elliptic system for the unknown (u%, u 2 )- In two dimen- 
sions, this system turns out to have a variational formulation with a coercive bilinear form: 

V • (A 2 \H\^H ij V Uj ) = 0, Ul \ dx = 9i . i = 1,2, (18) 

It thus admits a unique solution in the following space (up to an additive i^-lifting of {91,92)) 

W := (Hq(X)) 2 , u=( Ul ,u 2 )£n iff ||u|&:=/ \V Ul \ 2 + \Vu 2 \ 2 dx < 00. (19) 

Using a Fredholm-type argument, we obtain that this solution is also stable with respect to 
changes in the data, as stated in the following result: 

Proposition 2.5 (Stability of the strongly coupled elliptic system). Let H,H' have their com- 
ponents in W 1,oa (X) and satisfy ([7]). // u, u' are the unique solutions to (|18p with the same 
illumination g, then u — u' G H and satisfies the stability estimate: 

Hu-u'Hk < C\\H-H'\\ w ±,oo. (20) 

Once the couple (111,112) is reconstructed throughout X, one may reconstruct l^l" 1 using (a 
modified version of) the gradient equation ([8]). 

Theorem 2.6 (Stability of |^4|, second approach). Let the conditions of proposition \2.5\ be 
satisfied. Then the reconstructed determinants \A\ and \A'\ satisfy the stability estimate: 

\WA\~ 1 - lA'r^HHx) < C\\H - H'\\ w ^ {x) . (21) 
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The rest of the paper is structured as follows. We first derive equations ([8]) and (jlOp in 
section [31 which form the cornerstone of all subsequent derivations. Section [J] covers the three 
reconstruction algorithms mentioned above. Section [5] provides a proof of lemma 12.11 while 
section [6] concludes with numerical examples for each reconstruction algorithm. 



3 Proof of equations [8] and ( fTOjl 
3.1 Geometrical setting and preliminaries 

In this section, we work on an open connected subset f2 Q X, over which (Si, £2) satisfy 
the positivity condition flTJ). For a matrix M = PDP T E M%{2) with P E 0(2,M),D = 
diag (Ai,A2), and a scalar r 6 I, we can define uniquely M r := PD r P T E M s K r(2) by taking 
the positive r-root of each eigenvalue. Now, because A r = 72 is uniformly elliptic for any r£R, 
the vector fields (A r Si, A r S 2 ) also form an oriented frame (denoted A r S). The measurements 
can also be represented as 

H ij =A r S i 'A- r Sj, l<i,j<2, rGl. (22) 

From this assumption, one can deduce that any vector field V can be represented by means of 
the formula 



V = H ij (V ■ A r Si)A- r Sj, r£l. 



(23) 



Both equations (I22p and (|23l) only "see" A up to a scalar multiplicative constant, thus these 
equations still hold if we replace A by A = \A\~z A. 

Finally, we give the following important relation (only valid when n = 2) true for any 
M E M S {2,M): 



MJM = (detM)J, J :-- 



-1 

1 



(24) 



3.2 Proof of the gradient equation (jHJ) 

The derivation relies on the analysis of the properties of the vector fields J A" 1 Si for i = 1,2. 
First notice that since J is skew-symmetric, we have 

JA^Si- A~ 1 S i = 0, 1 = 1,2, 
Then, using the relation fl2U) with M = ^l -1 and the fact that J Si ■ S 2 = det(Si, S 2 ) =: , 
JA" 1 ^ • A~ 1 S 2 = -JA~ 1 S 2 ■ A' 1 Si = (A- 1 JA^Si) ■ S 2 = J5i • 5 2 = lAl'^H^. 
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This means that the vector fields J A x Si can be expressed using the representation (f23|) with 
r = —1: 

JA^Sx = H pq (JA~ 1 Si ■ A' 1 S p )AS q = H 2q \A\~ l \H\^ AS q , 
JA~ l S 2 = H pq {JA~ 1 S 2 ■ A~ 1 S p )AS q = —H lq \A\~ 1 \H\^ AS q . 

We now apply the divergence operator to ([25]) . Together with the fact that V • (Jj4 _1 5j) = 
-(JV) • (A^Si) = and equation (ji]), and using the identity V(fV) = V/ • V + /V • V, we 
derive 

Vj^r 1 • \H\^H qp AS p + \A\~ 1 V(\H\^H qp ) ■ AS P = 0, q= 1,2. 
Multiplying the last equation by A~ 1 S q , summing over q and dividing by |^4| _1 |i/|2 , we obtain 

H qp (-V log \A\ ■ AS p )A~ 1 S q + |.ff|-3(V(|iZ'|*.ff 9p ) • AS p )A~ 1 S q = 0. 

The first term is of the form (|23p with r = 1 and V = — Vlog|^4| and thus equals — Vlog|A|. 
We obtain the second term of the r.h.s. of (|8|). The first term of the r.h.s. of (jHj) is obtained 
from the second one by expanding the term V(|iJ | 5 H pq ) and using the product rule and identity 
(f23|) to obtain the ^Vlog \H\ term. 

3.3 Proof of (HO]) 

We now orthonormalize S into a frame R via the transformation R = ST T , also written as 

R{ = t^jSji Sj = t ^ Rj , 1 ^ 2 ^ 72. 

The matrix T satisfies T T T = H -1 , also written as tkitkj = H l i for 1 < i,j < n = 2. It can be 
constructed by the Gram-Schmidt procedure or by writing T = H~2. With the V^-'s defined in 
@, the following important identity holds 

(V//'')/' 7 '/' 7 = (V(tpit pj ))t ik P l = V(V^ j7 + V(VV)^" = ^ + Vik, 1 <l,k <2. (26) 

Therefore, starting from (|8|), we have 

V log \A\=N+ (VH jl ■ ASi)A~ x Sj = N+ ((V H jl )t lp t jq ■ AR p )A~ 1 R q 

= N + ((V pq + V qp )-AR p )A~ 1 R q , (27) 

where we have used (j26|) in the last equality. Now, in order to derive equation (|10|) . we must 
write the Lie bracket [AR 2 , ARi] in two different manners. 
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On the one hand, writing [AR2, ARi] in the canonical basis (ei,e2) and using the identity 
[aX, bY] = a(X ■ V)(b)Y - b(Y ■ V)(a)X + ab[X, Y], we have that 

[AR 2 , ARi] = [AijRie^AkiRiek] 

' A ij A kl (R j 2 d i R l 1 - R{diR 2 ) + AijdiAuiRiRi - R{R l 2 )) e k , 



after renumbering indices properly. Moreover, in the parameterization R(0) we have 

r * q - f ■ = 7 f if 3 = I 

R{m\-R[d % R l 2 = \ 7 I/" and R{R\ - R[R l 2 = { 1 if (j, Z) = (2, 1) , 
^ I -1 if (J,Z) = (1,2) 

thus we obtain 

[AR 2 , ARi] = ((AnAki + A i2 A k2 )di9 + A i2 ^2 fc i - A^A^) e k = A 2 V8 + [A 2 ,A X \. (28) 



On the other hand, we compute [AR 2 , AR\] using @. First, deriving a divergence equation 
for ARi, i = 1, 2, and using the fact that ([U) can be recast as V • (ASi) = — log \A\ • ASi, we 
have 

V • (ARi) = V • (AtijSj) = Vtij ■ ASj + %V • (AS^) 

= • AV k R k - Uj^V log \A\ ■ ASj = V lk ■ AR k - log \A\ ■ ARi 

= A N .AR l + Vg-AR k , (29) 



where we have used (J27J) in the last equality. We now use the following 2d vector calculus identity 
[U, V] = (U ■ V)V - (V • V)U = V ■ (V <g> U - U <g> V) - (V • U)V + (V • V)U, (30) 
where V • M := diM^e^ for M a 2 x 2 matrix. With U = AR 2 and V = AR\, we have 

ARi ® AR 2 - AR 2 <g> ARi = A(R 1 ®R 2 -R 2 ® Ri)A = -AJA = -(det A) J = -J, 
so the first term in the r.h.s. of (1301) is zero. Thus we have 



[AR 2 ,AR 1 ] = (V • ARJAR2 - (V • AR 2 )ARi 

= --(N ■ ARi)AR 2 + -(N ■ AR 2 )ARi + (Vfe ■ AR 2 )AR 2 + (V& ■ AR 1 )AR l 
= A(Rx ®R! + R 2 ® R 2 )AV 1 a 2 - ^A(R 2 ® i?i - R x ® # 2 )liV 

where we have used the properties R\ ® R\ + R 2 ® R 2 = I 2 and R 2 £3 i?i — i?i <8> R2 = 
Combining fl28f) with the last r.h.s. yields (I10p . 



4 Reconstruction procedures 

This section is devoted to the reconstruction procedures. We first reconstruct the anisotropic 
part of the conductivity tensor A and second reconstruct (#,log \A\) and {u\,U2, \A\~ 1 ). 

4.1 Reconstruction of the anisotropy A = 72 with m = 3 or 4 

Let us now consider a quadruple (<?i, <? 2 , 53, #4) £ G-y with possibly g 2 = #3- Condition (JTTD 
ensures that the matrices = [s[ 1] \S^] = [Si\S 2 ] and = [S[ 2) \S^ 2) ] = [S 3 \S A ] satisfy the 
positivity condition (|7|). Let us denote i?W = S®TW T the SO(2, R)-valued matrix obtained 
after Gram-Schmidt orthonormalization, parameterized by an angle function 0j, and denote 

H (i) =3 (iyr s (i) i N (i) : =Iviog|fl-»|, : = _ 4 = 1,2. 

For each pair of solutions, we have the equation 

1 2 W 4 + [1 2 , lx] = &vf® - Y - JJVW, i = 1,2. (31) 

Taking the difference of both systems cancels the term L42,Ai], and we obtain equation (|14p . 
with vector fields 

= (*i,* 2 ) T := V(fl 2 - Oi) - V$ 2) + and Y g = ( yi ,y 2 ) T := -\j{N^ - iV«). 

Now we claim that although the angle functions 61,62 are unknown, V(# 2 — 6\) is known from 
the data. Indeed, we have that 

V(# 2 - 0i) = cos(0 2 - #i)V(sin(# 2 - 00) - sin(# 2 - 0i)V(cos(0 2 - 6q)), 

and then, 

cos(0 2 - 00 = R? • flf> = t^S? ■ Sf\ sin(0 2 - e x ) = f#> ■ R? = t^s? ■ sf\ 

where both r.h.s. only depend on the data. As a result, the vector fields X s , Yg are completely 
known from the data {Hij}jj =1 . 

Parameterization of A 2 and inversion: The matrix A 2 is symmetric and has unit deter- 
minant and as such is characterized by only two scalar parameters. This is why equation (|14[) 
is enough to reconstruct it algebraically wherever X s 7^ or Yg 7^ 0. We now parameterize A 2 
as follows 



A 2 (£X) 



C 



1+c 2 



£ > 0, (32) 
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where the second row is deduced from the first one by constructing a symmetric matrix with 



unit determinant. With X, 



g 



(,Xg, Xg 



2 \T 



and Yg = (y^,y"i) T , equation (fT4"j) becomes 



ygj wg; 

and + (1 + C 2 )^ = ^ 



which can be rewritten as 



y\ 



and can therefore be inverted as 

^(Xg.Yg^^ + tXg) 2 ) 



and C = (X s • 



1 W 



x g x g 



) 



(33) 



(34) 



Proof of theorem 



The only important point is to show that X s ■ Yg is never zero. Since 
condition (|12[) is satisfied, then Yg never vanishes over X. Since A E we have X s ■ Y s = 

WA^YgW 2 > Ko 2 ||Ygll 2 - Therefore, X s ■ Y g = wher ever Yg = 0, that is, nowhere in this case. 
Inequality (|16p follows provided that the inequalities (|15|) hold in the Gram-Schmidt case, and 
the expressions (|34p are smooth in the components of X S ,Y S away from {X s ■ Yg = 0}. □ 

Remark 4.1 (An unstable parameterization of A). Another way (seemingly geometrically more 
meaningful) to parameterize A is, for (a, a) G (0,oo) x S 1 , to write 



A(a, a) 



c a s a 




' a 






s a c a 




a' 1 




Sq, Cq, 



(cq,Sq) := (cos a, since). 



(35) 



a describes the anisotropy and a locates the main axes of the ellipse. However, besides the fact 
that this representation is not injective (indeed we have A(a, a) = A(a, a+ir) = A(a~ 1 , a+ir/2) ), 
the reconstruction of a becomes unstable as a approaches 1 . 



4.2 Reconstruction of \A\ 



i 1 1 

T 2 



We now consider the problem of reconstructing \A\ from m = n = 2 measurements, assuming 
that A is known. We assume that the positivity condition is satisfied throughout the domain 
X. We propose two approaches, which we analyse consecutively in the next two sections. 



4.2.1 Reconstruction of (6, log \A\) 

This approach consists in reconstructing 9 and then \A\ using the gradient equations (|10p and 
©. We first isolate V8 in (fTUj) by writing 



W = V& - A 



-JN+[A 2 ,A!] 



(36) 
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We thus require an expression for [y^j-^i]- I n the case where A 2 was reconstructed in the 
previous section using the (£, C) variables, we need to compute A from knowledge of A 2 , which 
we do by first introducing a parameterization of A similar to (|32p . called A(X, /j,) with A > 0. 
Then, equating the terms in the first row of both representations (A(X, /j,)) 2 and j4 2 (£,£), we 
obtain the relations 

X 2 + n 2 = £ and ^(l + A 2 + /i 2 ) = C, 
A 

which, using the condition A > we invert as, 

A = (1 + 6 (c 2 + (i + 6 2 ) 2 and ^ = C (c 2 + (i + 2 

In the variables (A,//), we now obtain the following expression for the term [^42,-Ai] 

i+M 2 



(37) 



[A 2 ,A 1 



— P — 



VA 



A 



/' 



u ^ 

r A 



V/x. 



(38) 



Using (|38p in (|36p allows us to reconstruct 9 via integrating the gradient equation (|36p along 
curves (and assuming that 6 is known at one point of the domain). Alternatively, if is known 
on the whole boundary, one may apply the divergence operator on both sides of (|36p and solve 
a Poisson equation with Dirichlet boundary conditions. 

On to the reconstruction of \A\, we now use equation ([8]) in the R(0) frame to obtain: 



1 2 _ „ 

V log \A\ = -V log \H\ + 2 ]T (t£ • ARp)A -1 i? g , 



(39) 



p,g=l 



whose r.h.s. is completely known. For its resolution, equation (|39p may be simplified using a 
calculation similar to [Ml Sec. 6.2]. To this end, we define &ij(6) = Ri® Rj for 1 < i, j < 2 and 
compute 



1 2 „ 

Vlog|A| = -Vlog|#| +2 A~ l % q AV; q 



p,q=l 

= -V n - V 22 + 2AT l <S> 11 AV 11 + 2A~ 1 <S> 22 AV 22 + A~ l {^ 12 + § 2l )A{V 12 + V 21 ) 
= I" 1 ^!! - *2a)A(Vii - V 22 ) + A7\* n + § 2 i)A(y 12 + V21), 

where we have used the facts that A~ 1 (&u + & 22 )A = I 2 and — V11 — V 22 = Vlog|i?|2. The 
matrices <3>ii — $22 and &12 + ^21 axe symmetric matrices that can be expressed in the following 
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manner 

$11 - $22 = cU + sJU and $12 + $21 = -sU + cJU, where 
(c,a) := (cos(20),sin(20)), U := 



1 
-1 



(40) 



From this we deduce the final expression of Vlog \ A\: 

V log \A\ = A- 1 ( cos(2fl)F c + sin(20) JF C ) , F c := UA(V U - V 22 ) + JUA(V U + Vn). (41) 

Proof of Theorem \2.4\ The proof is very similar to Theorem 3.2]. For two sets of measure- 
ments H and H' , the error function on 9 satisfies 

v(0 - e') = v? 2 - - \a~ 2 j{n - n'). 

Assuming that 6(xq) = 9'{xq) for some xq £ X and using Gronwall's lemma along (bounded) 
paths joining any x G X to Xq, and provided that ||Vij — V^||z/x> < C\\H — H'\\wi,oo in the 
Gram-Schmidt case, we arrive at the inequality 

\\9 - 9\\ w i, X ( X ) < C\\H - H \\ w i,ao( X )- (42) 
Similarly, using the difference of equation (|39p for log |^4| and log \A\' and using ()42p . we arrive 

at dm). □ 



Remark 4.2. As previously pointed out in | jgj/ . solving gradient equations may require the 
enforcement of compatibility conditions (i.e. the r.h.s must be curl-free), in order to ensure that 
the computed solution does not depend on the choice of integration curve. 

4.2.2 Reconstruction of (ui,U2,\A\~ 1 ) 



This approach, first introduced in [16], consists in writing a strongly coupled elliptic system 
for the unknowns (7/1,7/2), whose properties are particularly appealing in two dimensions. We 
proceed as follows. 

In the frame (Viti, V7/2), equation (jHJ) reads 

Vlog |A| = \A\\H\-5(V(\H\3H p *) -A 2 Vu p )Vu q . (43) 

Now, the diffusion equation ([3]) can be rewritten as 

V • {A 2 Vui) + V log \A\ • A 2 Vui = 0. 
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Plugging ([35]) into the latter equation and using the fact that l^jVuq • A 2 Vui = H q i yields the 
system 

V • (! 2 Vui) + W ip ■ A 2 Vu p = 0, Ui\ dx =g u i = 1, 2, where (44) 
W ip :=H gi \H\-^V(\H\^H^), l<i,p<2. (45) 

Upon multiplying (I44p by l-f^ai?- 7 ' and summing over j, we obtain an equivalent formulation in 
divergence form 

V ■{\H\hw i A 2 Vu i ) = 0, u i | ex = ^, j = 1,2. (46) 

We assume that 31,32 £ H?{dX) and define u$ to be a i^-lifting of inside X. Defining the 
new unknown w = (w±,W2) ■= (u\ — V\,U2 — V2), w satisfies the following two equivalent systems 
whenever u = (u\,U2) satisfies (1331) or ([36]) and vice-versa 

V • (A 2 Vwi) + W ip ■ A 2 Vw p = hi := -V • (A 2 V Vi ) + W ip ■ A 2 Vv p , Wi \ dx = 0, i = l,2. 

(47) 

-V-(\H\ 1 2A 2 H ki Vw i ) = f k :=V-(\H\^A 2 H kl Vv l ), x e X, = 0, A: = 1,2. 

(48) 

System (|48p allows us to establish existence and uniqueness of w while (|47p is used to establish 
the stability of w with respect to the data H . 

Uniqueness of (u\,U2)' System (|48p can be recast as the variational problem of finding 
w € % := Hq(X) 2 such that for every w' G we have 

B(w,w')= f f k w' k dx, where B(w,w') := f \H\ ^H ki (A 2 V Wi ) ■ Vw' k dx. (49) 
Jx J X 

With % endowed with the norm ||w||^ = f x \Vwi\ 2 + |Vw2| 2 dx, assuming the positivity con- 
dition ([7]), the matrices H, H^ 1 and A 2 are all uniformly elliptic over X, which guarantees that 
the bilinear form is coercive. The continuity of B and of the linear form w' 1— > j x f k w' k dx over 
T~L are straightforward. As a result, Lax-Milgram's theorem establishes existence and uniqueness 
ofweH solving (gSJ and ([37| . and thus ofuGv + 'H solving ([33]) and ([35]) . 

Stability of (ui,U2) with respect to the data: Let us define the operator L _1 : L 2 (X) B 
f 1— >■ w, where w is the unique solution to the problem 

Lw ■= -V • (A 2 Vw) = f in X, w\ dx = 0. (50) 
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Since A 2 is uniformly elliptic, the operator L _1 : L 2 (X) — > H 2 (X) is bounded (see e.g. [TT]). 
thus compact L 2 (X) — > Hq(X) by the Rellich compactness theorem. Therefore, applying 
to (|47p . we obtain the integro-differential system 

Wi + L'^-Wip ■ A 2 Vw p ) = -L~ x hi, i = 1,2, 

which in vector notation may be recast as 

(I + P w )w = {L~ 1 f l } 2 =1 , where P^/w = W ip • l 2 ^)}^ . (51) 

Similarly to [16l Lemma 5.1], the operator : % — > 7i is compact and its operator norm 
satisfies 

||PwHI < C||Wl|oo, Halloo = max WWijWwx), (52) 

l<i,j<2 

see [16] for details. Therefore, equation (|5ip is a Fredholm equation and the boundedness of 
(J + Pu--) -1 holds as soon as —1 is not an eigenvalue of Pw- This is the case here, since the 
previous paragraph proves exactly this fact. On to the stability of u w.r.t. the data H, we use 
the fact that the vector fields W defined in (j45|) satisfy estimates of the form 

l|W - W\\ L oo( X ) < C\\H - H\\ w i,oo( X )i 

whenever H and H' have their components in W l >°°(X) and satisfy the positivity condition ([7]). 
With the previous estimate, the proof of proposition 12.51 is similar to Proposition 2.6] so we 
do not reproduce it here. 

Reconstruction of | ~ 1 : On to the reconstruction of \A\, equation (I43p can be recast as 

V\A\~ l = -\H\-^(V(\H\tH pq ) ■ A 2 Vu p )Vug, (53) 

the r.h.s. of which is now completely known. One may thus choose to solve this equation either 
solving ODE's along curves or taking divergence on both sides and solving a Poisson equation. 
The stability of such a reconstruction scheme was already stated in theorem 12.61 whose proof 
may be found in the very similar (isotropic) context of [W\ Theorem 2.8]. 

5 Proof of lemma 12.11 

Isotropic case 7 = CJI2: Consider the problem V • a(x)Vu = on 1R 2 with cr(x) extended 
in a continuous manner outside of X and such that a equals 1 outside of a large ball. Let 
q{x) = on M 2 . We assume that q G H 3+£ (M. 2 ), which holds if a - 1 € H 5+e (R 2 ) for some 
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e > 0, i.e., the original a\ x £ H 5+£ (X). Note that by Sobolev imbedding, a is of class C 4 (X) 
while q is of class C 2 (X). 

Let v = y/au so that Av + qv = on M. 2 . Let p G C 2 be of the form p = p(k + ik -1- ) with 
k G S 1 and k- 1 = Jk so that k • k- 1 = 0, and p = \p\/y/2 > 0. Thus, p satisfies p ■ p = and 
e p x is a harmonic complex plane wave. Now, it is shown in [8], following works in [91 [18], that 

u p = VoTip = e"- !B (l + Vp), (54) 

with (A + q)v p = and hence V • aX7u p = in M 2 . Furthermore, with the assumed regularity, 
[HJ Proposition 3.3] shows the existence of a constant C such that 

P\\^p\\c 2 (x) ^ Clkllfra+spQ, an d so lim^ ||^p|| C 2 ( x) = °- (55) 
Taking gradients of the previous equation and rearranging terms, we obtain that 

V^Vup = e p - x (p + ip p ), with ip p := V?p p + ip p p - (1 + t/; p )Vy/a. (56) 

Note that u p is complex-valued and since a is real-valued, both the real and imaginary parts 
and «p are solutions of V • (uVn) = 0. More precisely, we have 

v^Vn^ = pe^ x ( (k + p" Vp) cos(pk ± • x) - (k x + p~Vp) shi^k 1 • xj) , 

(57) 

v^Vn^ = pe pkx ((k x + p~V?) cosCpk 1 • x) + (k + p~ Vp) single 1 • x)) . 

We now denote d p := det(y/aVu^, ^/aVu^), and straightforward computations lead to 

d p = pV**(l + / P ), / p := p- 1 (k-L • <pj + k • *,*) + P ~ 2 J^ P ■ <p% 

where lim^oo sup^ |/p| = 0. Letting p so large that sup^l/pl < |j the function d p is now 
bounded away from zero over X. Taking logarithm and gradient, we now obtain 

Vlogd^pfWp^ 1 ^--) . (58) 

We now define ki,k2 £ S 1 such that ki ^ k2 and define, for j = 1,2, pj := p(kj + ikj~). 
Considering the solutions {u pi , tt^, w^ 2 , itp 2 ), the previous calculations show that, for p large 
enough, we have 

inf min(d Pl , d P2 ) > Co > 0, 
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which means that condition (jlip is satisfied. Furthermore, using (|58p . we have 



d Pi _ „mn, -1 / V ^ v /w 



Vlog-^=p(2(k 1 -k 2 )+r), r:=p 



dp 2 V 1 + fpi 1 + /pa 



where the remainder r may be made negligible by virtue of (|55|) and the smoothness assumption 
on a. For p such that sup^ ||r|| < ||ki — k 2 ||, the quantity V \og(d pi /d p2 ) never vanishes over X 
and thus condition (j!2p is satisfied. In this case, let g p = {c/i}i<i<4 be the traces of the above 
CGO solutions (it^ , tt^ , iip 2 , ) . These illuminations generate solutions that fulfill conditions 
(fTTI) and (fT2j) , so g p G ^ CT i 2 . By continuity, any g in an open set (of sufficiently smooth boundary 
conditions) in the vicinity of g p also belongs to Q a \ 2 and the isotropic case is proved. 

General case: Since 7 = |7| _ 27 is a K 2 -conformal structure on C, Theorem 10.2.2] implies 
that there exists a unique diffeomorphism (f> : C 3 z 1— > </>(z) = z' € C satisfying the Beltrami 
system (normalized at infinity) 

D t (t}jo^D(/) = J (t> h, := det(Ity), z G C, 0(z) z =°° z + C^z" 1 ), 

alternatively formulated as the following complex Beltrami equation 

d $ (a.i \\ d ^ v, 722 - 7n- 2 *7i2 - 
= u {^ z )) where v= — , z G C. 

oz dz 2 + 711 + 722 

v defined above coincides with (JT3J) and thus belongs to Cf oc by assumption. By virtue of [21 
Theorem 15.0.7], this implies that is a C^-diffeomorphism. We further have Ja, > c\ > 
throughout C. Using this change of variables and denoting V = ^7 x ,y an d V = V^/y in the 
sequel with x' + iy' = z' and x + iy = z, it is well-known that a function u solves 

V • ( 7 Vu) = 0, zeC, (59) 

if and only if the function v = u o 4>~ l solves 

V • (<P^Vv) = 0, where ^ 7 (z') = — ^L><^(z) 7 (<^)) £><K*) = °(z')h, (60) 

with cr(z') := |'y] 2 o (/> -1 (z'). Using the fact that I7I2 G ff 5+e (X) by assumption and G C? , 
the change of variable for Sobolev spaces implies that a G H 5+£ (cj)(X)) . Thus by virtue of the 
first part of the proof, Q a i 2 7^ 0. 

Let g G Gah defined on the boundary <p(dX) = d{(j)(X)). Then it is easy to see that 
the illumination g o (j) G G-y, so that 7^ 0. Indeed, if for 1 < i < 4, Vi designates the 
unique solution to ([60]) over 4>{X) with boundary condition Ui|g(^m) = 9i, then the function 
Ui := Vi o (p satisfies (|59p over X with boundary condition = 9i 0- Using the chain rule 
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Vu = V(v o (p) = Dcf)(X7'v) o (ft and properties of the determinant, routine computations yield 
the following relations, true for every z = x + iy G X 

det(Vui,V Uj )(z) = J (j> (z)det(y'v l ,Vv J )((ft(z)), e {(1,2), (3, 4)} 

JY so4> (z) = D<j>(z) JY g (</>(z)), 

with Yg defined in f)12[) . Since D</> is everywhere invertible with > cq > 0, the previous 
relations imply that (u\, U2, U3, 114) satisfy conditions (fTT]) - (fT2]) over X if and only if (i>i, V2, i>3, V4) 
satisfy (fnj) - (fT2j) over <ft(X), that is, g S Q a i 2 if and only if g o (ft e Q^. The lemma is proved. 

Remark 5.1 (On the regularity of a). The existence of CGO solutions can be established in two 
dimensions assuming mere boundedness of the isotropic diffusion coefficient a, as was established 
in J3y. However, due to the necessity of estimate (|55|) for our purposes, we need the results in 
J2J/, which in turn require higher regularity on a. 



6 Numerical examples 

In order to validate the reconstruction algorithms presented in the previous sections, we imple- 
mented a forward solver for the anisotropic diffusion equation on a two-dimensional square grid, 
using a finite difference method implemented with the software MatLab. 

We use a N+lxN+1 square grid with N=128, the tensor product of the equi-spaced subdi- 
vision x = -l:h:l with h = 2/N. Partial derivatives are performed using the operators DX = 
kron(D,I) for d/dx and DY = kron(I,D) for d/dy, where D designates the one-dimensional fi- 
nite difference derivative matrix and I=speye(N+l) is the N+lxN+1 (sparsified) identity matrix. 
D is the second-order centered stencil [-1 l]/(2h) with, at the boundary D(l,l:3) = [-3 
4 -l]/(2h) and D(N+1 ,N-1 : N+l) = [1 -4 3]/(2h). 

In the following examples, the conductivity tensor is described by the triplet of scalar func- 
tions (|7|2,£,£) such that 7 = |tI ^ tC^; C) with 7 given in ([32]) . Note that \A\ = |-y 1 2 . The 



1 



anisotropy coefficients (£, C) in Fig. Wajfc (b) are used in all experiments, while the determinant 



I71 2 may be either one of the two functions displayed in Figs. Q ||c)p L; (d) 

We sometimes perturb the internal functionals Hij with random noise of the form 

OL 

H n oisy = H .* (1 + — random), (61) 

where a is the noise level and random is a N+lxN+1 matrix of random entries with uniform 
density over [—1,1], to which we have applied a slight regularization by convolving it with the 
averaging filter ones (3) /9 (e.g. using the MatLab imfilter function). 

Reconstruction of the anisotropy A from noiseless data and m = 3 solutions. Let 

the two systems = (Si,^) and = (82,83) be associated with (91,52) and (92,93) 
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(a) £ (b) £ (c) |7|a smooth (d) [-/I a rough 

Figure 1: The coefficients used in the numerical simulations 



respectively, where (<?i, 52,53) are given by 

{9i,92,93)(x,y) = (x + y,y + 0.1y 2 ,-x + y), (x,y) G d[-l,l} 2 . 
We first compute the data {Hij} by solving the three forward problems ([3]) and then computing 



Hi 



7Viij • Vuj over the grid. We then compute the vector fields X and Y in the Gram- 



Schmidt case, where the transfer matrices = {t fj} 2 j =1 such that = S^T^ T for I = 1, 2 
are given by 



T (i) 



and T( 2 ) 



n 22 







H 22^-2 



with di := (HnH 2 2 - Hf 2 )? and d 2 = (F 2 2#33 - #23) 
expressions: 



-H23H 22 2 d 2 1 
X g and Yg then admit the following 



H- 



22 



H 



12 



H; 



22 



4 V 



23 



H 



22 



and K 



g 



--JV(logd 2 -logdi). (62) 



Once Xg, Y s are generated numerically, we then reconstruct the functions (£, £) that characterize 
via formulas (1341 ). 

The simulation of log d\ — log d 2 that appears in (|62|) is presented for the smooth 7 in the 
absence of noise on Fig. E j](a)| and in the presence of one percent of noise on Fig. [ ^b)[ 

The reconstruction of (£, C) is performed for both forms of |-/| 2 in Figs{ ^c)^|(d)| and the 
results are presented in Figj2j In the smooth case, the reconstructed £, C cannot be distinguished 
visually from the exact coefficients and are thus not represented. Relative L 2 (L°°) errors for £ 
and ( are 0.1% (8.6%) and 0.8% (13.7%), respectively. In the second case, the singularities of 



I7I2 create artifacts on the reconstructions, see Fig. [ ^b)| fc (f) and the relative errors for £ and 
C increase to 5.4% (99%) and 15.8% (167%), respectively. 
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(a) true £ 



(e) true £ 



(bH("l» on(d)l 



(c) £ ("2" on (d)l 




(f) C ("1" on (h)l 



(g) C ("2" on (h)l 



— <L 2 
—*— % true 



,1 

j 'I; >* 1 



(d) £ at {x = 0.5} 




-1 -0.5 



0.5 1 



(h) C at {z = 0.5} 



Figure 2: Anisotropy reconstructions, (b)^(f) with rough j'yj a and noiseless data, (c)^(g) 



with smooth I71 2 ; noisy data (a = 0.1%) and p = 100 measurements. 



Reconstruction of the anisotropy A from noisy data and m 

Figs. 



3p solutions, p > 1. 



(b) show that even very small amounts of noise in the three available internal 



functionals significantly affect the reconstruction of the anisotropic coefficients. The presence 
of noise creates local extrema for the function logdi — logc^ so that its gradient and Y s may 
vanish and (|34h may no longer hold. 

To address this lack of robustness, we increase the number of available internal functionals 
to 3p for p > 1, which correspond to the boundary conditions (gi, . . . , g p ). Instead of solving 
the linear system (|33p . we solve the normal equation (in the L 2 -minimizing sense) associated 
with the over-determined 6p linear equations, each pair of which looks like (|33p with g = gj for 
1 < j < P- The normal system reads 



E 







E 

i=i 



«) 2 + (4) 2 
12 12 

X Si X Sj ~ y sj y s, 



sym 

i) 2 + «) 2 
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121 • 



(63) 



and thus may be inverted as 

p 

(£, C) = (det Z)- 1 (Y, X zr Y z) ( H 22, -H 2 
i=i 

Results are shown after p = 100 iterations in Fig. [2j where for 1 < j < p, we used the following 
boundary conditions 

giJ (x,y) = (3+(x,y)-P) 1+ i/P, /3 := := (cos A sin/3), 



7T 



g 2 j{x,y) = (x,y) • /3 + - + 0.01 2+(z,y) •£ + 



7T 



\ 93,j(x,y) = \ 3 + {x,y) ■ f3 + 



7T \ P 



The relative L 2 errors for £ and £ are 22% and 27%, respectively. The effect of the number p of 
measurements on the reconstruction can been on Fig. 02 We observe a very slow convergence, 
which is consistent with the central limit theory, as p increases. The slow convergence may 
be considerably sped up by adding more constraining prior information on (£, £) such as for 
instance regularity or sparsity constraints. We do not explore this aspect here. 




-0.01 -0.005 0.005 0.01 




-0.015 -0.01 -0.005 0.005 0.01 0.015 



(a) log di - log d 2 (a = 0%) (b) log di - log d 2 (a = 1%) (c) £ at {x = 0.5} 



Figure 3: (a) (b) influence of the noise on the function log d\ — log d.2- (c) li (d) 
of £ and £ using reconstruction formulas (j63H for a few values of p. 



(d) C at {x = 0.5} 

cross sections 



Reconstruction of |-y| 2 via (0, log I'yl 2 )- We now assume that the anisotropy A is known 
or reconstructed, and solve for 6 by first applying the divergence operator to (|36p and then for 
log I71 2 by applying the divergence operator to (|4ip . For both Poisson equations, we use exact 

data as boundary conditions. In the Gram-Schmidt case, 6 = AS\ and the vector fields {Vij} 
are given by (recall that d := {H11H22 — Hf 2 )^) 

V n = VlogiT^, V 12 = 0, V2x = -^V0|^), V2 2 = V log^cT 1 ), 
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where the data come from solutions {111,112) with boundary conditions {gi, g 2 ){x,y) = (x,y) for 
(x,y) £9(-l,l) 2 . 




5 6 7 



(a) Hn (a = 0%) 



(b) Hu (a = 30%) 



-1 -D.5 0.5 1 



(c) H 12 (a = 0%) 



-1 -0.5 0.5 1 



(d) H 12 (a = 0%) 



Figure 4: Examples of measurement data. 



The isotropic part modeled by |-y| 2 is given in Fig. H ^e)[ Fig. H] displays the corresponding 
internal functionals with and without noise. Fig. [^displays the reconstructed 9 and I7I 2 . These 
reconstructions are quite robust to noise when the anisotropy is known: the relative L 2 (L°°) 
errors are 0.06% (0.14%) for |7|t and 0.04% (0.4%) for 9 with noiseless data, and 3.2% (12.5%) 
for |— y I 2 and 14.2% (24.5%) for 9 with 30% noise. If the anisotropy A is first reconstructed from 
noisy data, this certainly has repercussions on the reconstructed (9, I7I 2), as can be seen in Fig. 



c)p z (g) In this case, the L (L°°) relative errors are 2.2% (9%) for | 'y 1 2 and 42.7% (60%) for 



Reconstruction of | y 1 2 via {111,112, \"f\ 2 ). We still assume A known with the coefficients 
(£, C) displayed in Fig. [2j The boundary conditions are the same as in the preceding example. 
The isotropic component |-y 1 2 is now given in Fig. [ ^e)| and corresponds to a non-smooth 
coefficient. In this context, we first reconstruct (^1,1*2) by solving system (|46p . after which we 
reconstruct | y | 2 by taking the divergence of (|53p . The relative L 2 (L°°) errors are 3.9e — 13% 
(6.8e - 13%) for m, 2.5e - 13% (6.8e - 13%) for u 2 and 13% (62%) for | 7 |5 in the case of 
noisefree data, and 0.2% (0.7%) for m, 0.1% (0.4%) for u 2 and 14% (71%) for in the case 
of data polluted with 30% noise. See Fig. [6] for the display of some reconstructions. Based on 
the numerical simulations that we have performed, this reconstruction method and the previous 
one are very comparable in terms of accuracy and robustness. 
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(a) true 9 



(b) 6 ("1" on [(d)]) 



(c) 6 ("1" on [(d)]) 



1 -0.5 0.5 



(d) 6 at {x = 0.5} 





(e) true | 7 |' 



(f) | 7 |t ("1" on (h)l 




(g) | 7 |3 ("2" on (h)l (h) | 7 |3 at {a; = 0.5} 



Figure 5: Reconstruction of (6, |-y 1 2 ) . (b)k(f) with true anisotropy and noisy data (a = 30%). 



c)fc(g) with noisy data (a = 0.1%) using reconstructed anisotropy from FigMc)& (g 
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(a) tii 



(e) true | 7 |: 



(b) mi at {x = 0.5} 



(c) U2 



(d) U2 at {x = 0.5} 



(f) | 7 |t ("1" on (h)l 




(g) | 7 |2 ("2" on (h)l 



(h) |7|3 at {x = 0.5} 



Figure 6: Reconstruction of (u±,U2, \"f \ 2 ) with true anisotropy and discontinuous | 7 | 2 . |(f)[ with 
noisefree data, (g) with noisy data (a = 30%). 
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7 Conclusions 



This paper presented an explicit reconstruction procedure for a diffusion tensor 7 = (7r/)i<«,j<2 
from power density internal functionals H^ = ^yVui ■ V-Uj for 1 < i,j < I with / < 4. 

Provided that four illuminations g\ are selected so that the qualitative properties (jlip and 
(|12p of the corresponding solutions Ui are verified, we obtained an explicit expression for the 
anisotropic part of 7 and showed in Theorem l2.2l that errors in the uniform norm of the anisotropy 
were controlled by errors in the uniform norm of derivatives of the functionals H^. 

Once the anisotropy is reconstructed, or is known by other means, we have presented two 
methods to reconstruct the determinant of 7. The first functional is based on first reconstructing 
the angle 6 between e x and 7VU1. The second method is based on solving a coupled system of 
elliptic equations for (u±, U2). In both cases, we need that the three internal functionals Hn, H22 
and H\2 satisfy ([7]) in X. And in both cases, we obtain that the error in the uniform norm of 
the derivative of the determinant was controlled by errors in the uniform norm of derivatives of 
the functionals H^. 

This shows that the reconstruction of the determinant of 7 is more stable than that of the 
anisotropy of 7. Such a statement was verified by numerical simulations. In the presence of very 
limited noise generated by the numerical discretization, we obtained accurate reconstructions of 
the full tensor 7. However, even in the presence of quite small additional noise on the functionals 

, we observed that the reconstruction of the anisotropy was degrading very rapidly. On the 
other hand, the reconstruction of the determinant of 7, assuming the anisotropy known, proved 
very stable even in the presence of significant noise in the available functionals. In practice, this 
shows that appropriate regularization procedures on the anisotropy need to be introduced, for 
instance based on regularity or sparsity assumptions. This now standard step was not considered 
here. 

The functionals emerging from ultrasound modulation experiments are thus sufficiently rich 
to provide unique reconstructions of arbitrary diffusion tensors. This should be contrasted 
with reconstruction procedures based on boundary measurements of elliptic solutions, in which 
diffusion tensors are defined up to an arbitrary change of variables inside the domain of interest 
[3]. Moreover, reconstructions are stable with a resolution that is essentially independent of the 
location of the point inside the domain of interest. 

The reconstruction procedure presented here is two dimensional. Although this is not pre- 
sented here, the reconstruction of the determinant of 7 knowing its anisotropy generalizes to 
the ra-dimensional setting using techniques developed in [16] . The reconstruction of the full 
diffusion tensor in dimension n > 3 remains open at present. 
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